#install.packages("sf")
#install.packages("sp")
#install.packages("spdep")
#install.packages("spData") 
#install.packages("spatialreg") 

#install.packages("mapview")
#install.packages("ggplot2")
#install.packages("tmap")

Exploración de datos

Primero usar la libreria de sf para mapear y explorar nuestros datos

muni = st_read("./data/mx_city_metro.shp")
## Reading layer `mx_city_metro' from data source 
##   `/Users/irenefarah/Downloads/sp_econ/data/mx_city_metro.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 76 features and 40 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 2745632 ymin: 774927.1 xmax: 2855437 ymax: 899488.5
## proj4string:   +proj=lcc +lat_1=17.5 +lat_2=29.5 +lat_0=12 +lon_0=-102 +x_0=2500000 +y_0=0 +ellps=GRS80 +units=m +no_defs

Pueden usar attach(muni) para referirse a las variables por nombre, pero no lo hare en este caso.

Formato puede ser punto, línea, polígono.

Datos provienen de:

Nombre Tipo de variable Descripción
ID integer ID = estado + municpality código
ID_TEXT string ID (texto) = estado + municpality código
state string estado código
state_name string Nombre of estado
mun string municipio código
mun_name string Nombre of municipio
pop_10 integer Population in 2010
pop_15 integer Population in 2015
pop_20 integer Population in 2020
ppov_10 float Porcentaje de la población en pobreza en 2010
ppov_15 float Porcentaje de la población en pobreza en 2015
ppov_20 float Porcentaje de la población en pobreza en 2020
pepov_10 float Porcentaje de la población en pobreza extrema en 2010
pepov_15 float Porcentaje de la población en pobreza extrema en 2015
pepov_20 float Porcentaje de la población en pobreza extrema en 2020
peduc_10 float Porcentaje de la población con rezago educativo en 2010
peduc_15 float Porcentaje de la población con rezago educativo en 2015
peduc_20 float Porcentaje de la población con rezago educativo en 2020
pheal_10 float Porcentaje de la población con falta de acceso a servicios de salud en 2010
pheal_15 float Porcentaje de la población con falta de acceso a servicios de salud en 2015
pheal_20 float Porcentaje de la población con falta de acceso a servicios de salud en 2020
pss_10 float Porcentaje de la población con falta de acceso a seguridad social en 2010
pss_15 float Porcentaje de la población con falta de acceso a seguridad social en 2015
pss_20 float Porcentaje de la población con falta de acceso a seguridad social en 2020
pqdwel_10 float Porcentaje de la población con carencia de calidad a la vivienda en 2010 (calidad de pisos, techos, paredes y hacinamiento)
pqdwel_15 float Porcentaje de la población con carencia de calidad a la vivienda en 2015 (calidad de pisos, techos, paredes y hacinamiento)
pqdwel_20 float Porcentaje de la población con carencia de calidad a la vivienda en 2020 (calidad de pisos, techos, paredes y hacinamiento)
pserv_10 float Porcentaje de la población sin acceso a servicios básicos de la vivienda en 2010 (falta de acceso a agua potable, drenaje, electricidad, combustible para cocinar)
pserv_15 float Porcentaje de la población sin acceso a servicios básicos de la vivienda en 2015 (falta de acceso a agua potable, drenaje, electricidad, combustible para cocinar)
pserv_20 float Porcentaje de la población sin acceso a servicios básicos de la vivienda en 2020 (falta de acceso a agua potable, drenaje, electricidad, combustible para cocinar)
pfood_10 float Porcentaje de la población viviendo con inseguridad alimentaria en 2010
pfood_15 float Porcentaje de la población viviendo con inseguridad alimentaria en 2015
pfood_20 float Porcentaje de la población viviendo con inseguridad alimentaria en 2020
ping1_10 float Porcentaje de la población por debajo de la línea de pobreza en 2010
ping1_15 float Porcentaje de la población por debajo de la línea de pobreza en 2015
ping1_20 float Porcentaje de la población por debajo de la línea de pobreza en 2020
ping2_10 float Porcentaje de la población por debajo de la línea mínima de pobreza en 2010
ping2_15 float Porcentaje de la población por debajo de la línea mínima de pobreza en 2015
ping2_20 float Porcentaje de la población por debajo de la línea mínima de pobreza en 2020

Ver proyecciones de los datos:

st_crs(muni)
## Coordinate Reference System:
##   No user input
##   wkt:
## PROJCS["Lambert_Conformal_Conic",
##     GEOGCS["GCS_GRS 1980(IUGG, 1980)",
##         DATUM["unknown",
##             SPHEROID["GRS80",6378137,298.257222101]],
##         PRIMEM["Greenwich",0],
##         UNIT["Degree",0.017453292519943295]],
##     PROJECTION["Lambert_Conformal_Conic_2SP"],
##     PARAMETER["standard_parallel_1",17.5],
##     PARAMETER["standard_parallel_2",29.5],
##     PARAMETER["latitude_of_origin",12],
##     PARAMETER["central_meridian",-102],
##     PARAMETER["false_easting",2500000],
##     PARAMETER["false_northing",0],
##     UNIT["Meter",1]]
municipios_4326 = st_transform(muni, 4326)
st_crs(municipios_4326)
## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCS["WGS 84",
##     DATUM["WGS_1984",
##         SPHEROID["WGS 84",6378137,298.257223563,
##             AUTHORITY["EPSG","7030"]],
##         AUTHORITY["EPSG","6326"]],
##     PRIMEM["Greenwich",0,
##         AUTHORITY["EPSG","8901"]],
##     UNIT["degree",0.0174532925199433,
##         AUTHORITY["EPSG","9122"]],
##     AUTHORITY["EPSG","4326"]]

Mapas coropléticos

ggplot(municipios_4326, aes(fill = ppov_20)) + 
  geom_sf() +  # tells ggplot that geographic data are being plotted
  scale_fill_viridis_c() +
  theme_minimal() + 
  labs(title = "Porcentaje de pobreza, 2020")

# Mapa interactivo
mapview(municipios_4326, zcol = 'ppov_20')
# Agregar los dos
mapview(municipios_4326, zcol = 'ppov_20') +
  mapview(muni, color = 'red', alpha.regions = 0, legend=FALSE)

Guía de código de colores en R

#tmap_mode('plot') # estaticos
tmap_mode('view') # interactivo
## tmap mode set to interactive viewing
tm_shape(municipios_4326) + 
  tm_polygons(col = 'tan', 
              border.col = 'darkgreen', 
              alpha = 0.5) # transparencia

Examinar nuestros datos

tmap_mode('plot') # estaticos

pe = tm_shape(municipios_4326) + 
  tm_polygons(col = 'pepov_20', 
              border.col = 'darkgreen', 
              alpha = 0.5) + # transparencia
    tm_layout("Pobreza extrema, 2020",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)

ping1 = tm_shape(municipios_4326) + 
  tm_polygons(col = 'ping1_20', 
              border.col = 'darkgreen', 
              alpha = 0.5) + # transparencia
    tm_layout("Por debajo línea ingresos, 2020",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)

educ = tm_shape(municipios_4326) + 
  tm_polygons(col = 'peduc_20', 
              border.col = 'darkgreen', 
              alpha = 0.5) + # transparencia
    tm_layout("Educación, 2020",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)
  
salud = tm_shape(municipios_4326) + 
  tm_polygons(col = 'pheal_20', 
              border.col = 'darkgreen', 
              alpha = 0.5) + # transparencia
    tm_layout("Salud, 2020",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)

sbv = tm_shape(municipios_4326) + 
  tm_polygons(col = 'pqdwel_20', 
              border.col = 'darkgreen', 
              alpha = 0.5) + # transparencia
    tm_layout("Servicios básicos de la vivienda, 2020",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)

alim = tm_shape(municipios_4326) + 
  tm_polygons(col = 'pfood_20', 
              border.col = 'darkgreen', 
              alpha = 0.5) + # transparencia
    tm_layout("Inseguridad Alimentaria, 2020",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)

tmap_arrange(pe, ping1, educ, salud, sbv, alim, ncol=3, nrow=2)

Pesos espaciales en R

#Importar de archivo existente (abrir archivo)
#Pasar de "neighbor lists" a "weight lists" (listw object)
queen.nb_existente = read.gal("./data/mx_city_metro_q1.gal", region.id=muni$ID,override.id=TRUE) #spdep

summary(queen.nb_existente)
## Neighbour list object:
## Number of regions: 76 
## Number of nonzero links: 400 
## Percentage nonzero weights: 6.925208 
## Average number of links: 5.263158 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8  9 10 11 
##  2  3 11 13 16 12  8  5  2  3  1 
## 2 least connected regions:
## 15034 15061 with 1 link
## 1 most connected region:
## 15033 with 11 links
q1_existente <- nb2listw(queen.nb_existente)

class(queen.nb_existente) # Identificar el tipo de objeto
## [1] "nb"
class(q1_existente) # Identificar el tipo de objeto
## [1] "listw" "nb"
#k4.nb_existente= read.gwt2nb("./data/mx_city_metro_k4.gwt", region.id=muni$ID)
#summary(k4.nb_existente)
# queen weights:
queen.nb <- poly2nb(muni, row.names = muni$ID)
q1 <- nb2listw(queen.nb)
summary(queen.nb)
## Neighbour list object:
## Number of regions: 76 
## Number of nonzero links: 400 
## Percentage nonzero weights: 6.925208 
## Average number of links: 5.263158 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8  9 10 11 
##  2  3 11 13 16 12  8  5  2  3  1 
## 2 least connected regions:
## 36 59 with 1 link
## 1 most connected region:
## 35 with 11 links
summary(q1)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 76 
## Number of nonzero links: 400 
## Percentage nonzero weights: 6.925208 
## Average number of links: 5.263158 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8  9 10 11 
##  2  3 11 13 16 12  8  5  2  3  1 
## 2 least connected regions:
## 36 59 with 1 link
## 1 most connected region:
## 35 with 11 links
## 
## Weights style: W 
## Weights constants summary:
##    n   nn S0       S1       S2
## W 76 5776 76 32.83595 317.4801

nn = n*n S0 = suma global de los pesos

# rook weights:
rook.nb <- poly2nb(muni, row.names = muni$ID, queen=FALSE)
summary(rook.nb)
## Neighbour list object:
## Number of regions: 76 
## Number of nonzero links: 388 
## Percentage nonzero weights: 6.717452 
## Average number of links: 5.105263 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8 10 
##  2  4 10 17 15 10  7  7  4 
## 2 least connected regions:
## 36 59 with 1 link
## 4 most connected regions:
## 35 65 67 74 with 10 links
plot(muni$geometry)
plot(queen.nb, st_geometry(st_centroid(muni)), cex = 0.6, pch = 20, col = 'red', add = TRUE)
plot(rook.nb, st_geometry(st_centroid(muni)), cex = 0.6, pch = 20, col = 'blue', add = TRUE)

nc.5nn<-knearneigh(st_geometry(st_centroid(muni)), k=5, longlat=TRUE)
nc.5nn.nb<-knn2nb(nc.5nn)
summary(nc.5nn.nb)
## Neighbour list object:
## Number of regions: 76 
## Number of nonzero links: 380 
## Percentage nonzero weights: 6.578947 
## Average number of links: 5 
## Non-symmetric neighbours list
## Link number distribution:
## 
##  5 
## 76 
## 76 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 with 5 links
## 76 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 with 5 links
knn5 <- nb2listw(nc.5nn.nb)
plot(muni$geometry)
xy <- st_coordinates(st_centroid(muni))[,1:2]
plot(nc.5nn.nb, xy, pch=20, col="blue", add=T)

Autocorrelación espacial en R

La dependencia espacial puede expresarse mediante autocorrelación espacial. La autocorrelación espacial revela valores similares en ubicaciones cercanas; se puede estudiar global o localmente. Concretamente, la autocorrelación espacial revelará si existe aleatoriedad espacial en la distribución de la inseguridad alimentaria entre municipios.

Para estudiar distintos patrones espaciales, la definición de pesos espaciales es crucial. La estructura de pesos espaciales determina la conectividad entre ubicaciones vecinas.

Estadísticas de autocorrelación espacial global

Ejemplos: Indice de Moran, Geary C, Getis-Ord, entre otras (mezclar similitud de atributos con pesos espaciales). \[ \sum_{ij}= w_{ij}f(x_ix_j) \]

Variable espacialmente rezagada

#Examinar si hay autocorrelación espacial con variable de interés
muni$Vx <- lag.listw(q1, muni$pqdwel_20)
plot(muni$pqdwel_20, muni$Vx)

Moran’s I (la estadística más usada): Los mapas muestran patrones de distribución espacial, dando indicios de presencia de autocorrelación espacial. Para formalmente evaluar la presencia de autocorrelación espacial global, utilizamos el índice univariado global de Moran.

\[ I= [\sum_i\sum_j w_{ij} z_iz_j/S_0]/[\sum_iz_i^2/N] \] donde \(S_0=\sum_i \sum_j w_{ij}\) y \(z_i = y_i-m_x\) es la desviación del promedio.

Es una estadístics de producto cruzado \((z_i z_j)\) similar a un coeficiente de correlación. En este caso, los valores dependen de una estructura de pesos \((w_{ij})\)

Hipótesis nula: La distribución de los datos es aleatoria en el espacio.

Consecuentemente, cuando la hipótesis nula es rechazada, podemos decir que hay diferencias espaciales en nuestros datos.

La autocorrelación positiva indica clustering, la autocorrelación negativa indica dispersión espaial.

Hacer estudios con distintos tipos de pesos para evaluar qué tanto se sostienen nuestras pruebas.

Estadística que puede rechazar la hipótesis nula por encontrar autocorrelación, no-normalidad. No sabemos bien qué es lo que está mal con nuestros datos.

moran.test(muni$pqdwel_20, q1, alternative="two.sided")
## 
##  Moran I test under randomisation
## 
## data:  muni$pqdwel_20  
## weights: q1    
## 
## Moran I statistic standard deviate = 9.0233, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.644472659      -0.013333333       0.005314509
moranTest <- moran.test(muni$pqdwel_20, q1)
moranTest
## 
##  Moran I test under randomisation
## 
## data:  muni$pqdwel_20  
## weights: q1    
## 
## Moran I statistic standard deviate = 9.0233, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.644472659      -0.013333333       0.005314509
moranTest$statistic
## Moran I statistic standard deviate 
##                           9.023318
moranTest$p.value
## [1] 9.123843e-20

Estadísticas de autocorrelación espacial local

Los clústeres se pueden identificar a través de indicadores locales de asociación espacial (LISA, Anselin 1995). Estos indicadores muestran las desviaciones de los patrones globales generales al mostrar la contribución de cada observación a la autocorrelación global general, lo que permite una visualización de los puntos calientes/fríos de la variable de interés Como se mencionó anteriormente, diferentes ponderaciones estudian la presencia de patrones mediante el uso de distintas técnicas de conglomerados, como el I de Moran local univariado.

En esta sección vamos a identificar valores atípicos (outliers) mediante las estadísticas de autocorrelación espacial local, identificando hot spots y cold spots.

Los ejemplos anteriores de estadísticas de autocorrelación espacial global están diseñados para encontrar si los datos (en su totalidad) tienen patrones de aleatoriedad espacial. Sin embargo, no muestran la ubicación de los clústers.

\[ \sum_{j}=w_{ij} f(x_ix_j)\]

#Explicar cuadrantes
moran.plot(muni$pqdwel_20, q1)

local <- localmoran(muni$pqdwel_20, listw = q1)
#Mapear la estadistica local de Moran (Ii)
moran.mapa <- cbind(muni, local)
tmap_mode("view")
tm_shape(moran.mapa) +
  tm_borders(alpha=0.7)+
  tm_fill(col = "Ii",
          alpha = 0.5,
          style = "quantile",
          title = "Estadistica Local de Moran I")

Vecinos incluidos en el cluster.

#Examinemos aglomeraciones
quadrant <- vector(mode="numeric",length=nrow(local))

# Centra la variable de interés alrededor de promedio
m.pqdwel_20 <- muni$pqdwel_20 - mean(muni$pqdwel_20)     

# Centra estadística local de Moran's alrededor de promedio
m.local <- local[,1] - mean(local[,1])    

# Umbral de significancia
signif_1 <- 0.1 
signif_05 <- 0.05
signif_001 <- 0.001

# Construir cuadrantes para 0.1
quadrant[m.pqdwel_20 >0 & m.local>0] <- 4  
quadrant[m.pqdwel_20 <0 & m.local<0] <- 1      
quadrant[m.pqdwel_20 <0 & m.local>0] <- 2
quadrant[m.pqdwel_20 >0 & m.local<0] <- 3
quadrant[local[,5]>signif_1] <- 0   

# Armar mapa para 0.1
brks <- c(0,1,2,3,4)
colors <- c("white","blue",rgb(0,0,1,alpha=0.4),rgb(1,0,0,alpha=0.4),"red")
plot(muni$geometry,border="lightgray",col=colors[findInterval(quadrant,brks,all.inside=FALSE)])
box() #Agregar caja
legend("bottomleft", legend = c("No significativo","Bajo-bajo","Bajo-alto","Alto-bajo","Alto-alto"),
       title="p-value = 0.1", fill=colors,bty="n")

# Forma alternativa
muni$LISA=quadrant
colores <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusteres <- c("No significativo","Bajo-bajo","Bajo-alto","Alto-bajo","Alto-alto")
tmap_mode("plot")
p1 = tm_shape(muni) +
  tm_fill(col = "LISA", style = "cat", palette = colores[c(sort(unique(quadrant)))+1],
          labels = clusteres[c(sort(unique(quadrant)))+1], popup.vars = c("ID")) +
  tm_borders(alpha=0.5)+
  tm_layout("LISA pvalue = 0.1",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)
  

# Construir cuadrantes para 0.05
quadrant <- vector(mode="numeric",length=nrow(local))
quadrant[m.pqdwel_20 >0 & m.local>0] <- 4  
quadrant[m.pqdwel_20 <0 & m.local<0] <- 1      
quadrant[m.pqdwel_20 <0 & m.local>0] <- 2
quadrant[m.pqdwel_20 >0 & m.local<0] <- 3
quadrant[local[,5]>signif_05] <- 0  

muni$LISA=quadrant
p2 = tm_shape(muni) +
  tm_fill(col = "LISA", style = "cat", palette = colores[c(sort(unique(quadrant)))+1],
          labels = clusteres[c(sort(unique(quadrant)))+1], popup.vars = c("ID")) +
  tm_borders(alpha=0.5)+
  tm_layout("LISA pvalue = 0.05",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)

# Construir cuadrantes para 0.001
quadrant <- vector(mode="numeric",length=nrow(local))
quadrant[m.pqdwel_20 >0 & m.local>0] <- 4  
quadrant[m.pqdwel_20 <0 & m.local<0] <- 1      
quadrant[m.pqdwel_20 <0 & m.local>0] <- 2
quadrant[m.pqdwel_20 >0 & m.local<0] <- 3
quadrant[local[,5]>signif_001] <- 0  

muni$LISA=quadrant
p3 = tm_shape(muni) +
  tm_fill(col = "LISA", style = "cat", palette = colores[c(sort(unique(quadrant)))+1],
          labels = clusteres[c(sort(unique(quadrant)))+1], popup.vars = c("ID")) +
  tm_borders(alpha=0.5)+
  tm_layout("LISA pvalue = 0.001",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)

tmap_arrange(p1, p2, p3)

Con KNN = 5

Al considerar a los vecinos como todos los municipios abarcados dentro de una distancia determinada, los municipios más grandes tendrán menor conectividad y menos vecinos que los municipios más pequeños.

# con KNN = 5
local <- localmoran(muni$pqdwel_20, listw = knn5)

#Examinemos aglomeraciones
quadrant <- vector(mode="numeric",length=nrow(local))

# Centra la variable de interés alrededor de promedio
m.pqdwel_20 <- muni$pqdwel_20 - mean(muni$pqdwel_20)     

# Centra estadística local de Moran's alrededor de promedio
m.local <- local[,1] - mean(local[,1])    

# Umbral de significancia
signif_1 <- 0.1 
signif_05 <- 0.05
signif_001 <- 0.001

# Construir cuadrantes para 0.1
quadrant[m.pqdwel_20 >0 & m.local>0] <- 4  
quadrant[m.pqdwel_20 <0 & m.local<0] <- 1      
quadrant[m.pqdwel_20 <0 & m.local>0] <- 2
quadrant[m.pqdwel_20 >0 & m.local<0] <- 3
quadrant[local[,5]>signif_1] <- 0 

muni$LISA_kn=quadrant
colores <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusteres <- c("No significativo","Bajo-bajo","Bajo-alto","Alto-bajo","Alto-alto")
tmap_mode("plot")
p1 = tm_shape(muni) +
  tm_fill(col = "LISA_kn", style = "cat", palette = colores[c(sort(unique(quadrant)))+1],
          labels = clusteres[c(sort(unique(quadrant)))+1], popup.vars = c("ID")) +
  tm_borders(alpha=0.5)+
  tm_layout("LISA pvalue = 0.1",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)
  

# Construir cuadrantes para 0.05
quadrant <- vector(mode="numeric",length=nrow(local))
quadrant[m.pqdwel_20 >0 & m.local>0] <- 4  
quadrant[m.pqdwel_20 <0 & m.local<0] <- 1      
quadrant[m.pqdwel_20 <0 & m.local>0] <- 2
quadrant[m.pqdwel_20 >0 & m.local<0] <- 3
quadrant[local[,5]>signif_05] <- 0  

muni$LISA_kn=quadrant
p2 = tm_shape(muni) +
  tm_fill(col = "LISA_kn", style = "cat", palette = colores[c(sort(unique(quadrant)))+1],
          labels = clusteres[c(sort(unique(quadrant)))+1], popup.vars = c("ID")) +
  tm_borders(alpha=0.5)+
  tm_layout("LISA pvalue = 0.05",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)

# Construir cuadrantes para 0.001
quadrant <- vector(mode="numeric",length=nrow(local))
quadrant[m.pqdwel_20 >0 & m.local>0] <- 4  
quadrant[m.pqdwel_20 <0 & m.local<0] <- 1      
quadrant[m.pqdwel_20 <0 & m.local>0] <- 2
quadrant[m.pqdwel_20 >0 & m.local<0] <- 3
quadrant[local[,5]>signif_001] <- 0  

muni$LISA_kn=quadrant
p3 = tm_shape(muni) +
  tm_fill(col = "LISA_kn", style = "cat", palette = colores[c(sort(unique(quadrant)))+1],
          labels = clusteres[c(sort(unique(quadrant)))+1], popup.vars = c("ID")) +
  tm_borders(alpha=0.5)+
  tm_layout("LISA pvalue = 0.001",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)

tmap_arrange(p1, p2, p3)

Heterogeneidad espacial en R

Se espera heterogeneidad espacial en México, donde algunas regiones son más heterogéneas que otras debido a las diferencias geográficas y económicas.

Imponer estructura de forma discreta (regímenes) o continua (smoothing, GWR - regresión geográficamente ponderada).

#Crear una variable dicotómica que sea 1 si el municipio es parte de la Ciudad de México y 0 si no.
muni$cdmx= ifelse(muni$state=="09", 1, 0)
municipios_4326$cdmx= ifelse(municipios_4326$state=="09", 1, 0)

# Ejemplo de manipulación de datos
#municipios_4326%>%group_by(cdmx)%>%summarise(n=n())
ggplot(municipios_4326, aes(fill = cdmx)) + 
  geom_sf() +  # tells ggplot that geographic data are being plotted
  theme_minimal() + 
  labs(title = "Régimen espacial")

No sabemos por qué hay diferencias, pero resolvemos el problema de heterogeneidad.

#Crear una variable dicotómica que sea 1 si el municipio es parte de la Ciudad de México y 0 si no.
municipios_4326$dens=as.numeric(municipios_4326$pop_20/st_area(municipios_4326)*1000)

mean(municipios_4326$dens)
## [1] 4.209716
municipios_4326$pop_dens= ifelse(municipios_4326$dens>4.222332, 1, 0)

#municipios_4326%>%group_by(pop)%>%summarise(n=n())
ggplot(municipios_4326, aes(fill = pop_dens)) + 
  geom_sf() +  # tells ggplot that geographic data are being plotted
  theme_minimal() + 
  labs(title = "Régimen espacial - densidad poblacional")